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ABSTRACT 

We demonstrate the use of machine learning algorithms in combination with segmentation tech¬ 
niques in order to distinguish coronal holes and filaments in SDO/AIA EUV images of the Sun. 
Based on two coronal hole detection techniques (intensity-based thresholding, SPoCA), we pre¬ 
pared data sets of manually labeled coronal hole and filament channel regions present on the Sun 
during the time range 2011 - 2013. By mapping the extracted regions from EUV observations onto 
HMI line-of-sight magnetograms we also include their magnetic characteristics. We computed 
shape measures from the segmented binary maps as well as first order and second order texture 
statistics from the segmented regions in the EUV images and magnetograms. These attributes 
were used for data mining investigations to identify the most performant rule to differentiate be¬ 
tween coronal holes and filament channels. We applied several classifiers, namely Support Vector 
Machine, Linear Support Vector Machine, Decision Tree, and Random Forest and found that all 
classification rules achieve good results in general, with linear SVM providing the best perfor¬ 
mances (with a true skill statistic of « 0.90). Additional information from magnetic field data 
systematically improves the performance across all four classifiers for the SPoCA detection. Since 
the calculation is inexpensive in computing time, this approach is well suited for applications on 
real-time data. This study demonstrates how a machine learning approach may help improve upon 
an unsupervised feature extraction method. 

Key words. Solar wind - Coronal holes - Filament channels - Feature extraction - Supervised 
Classification - Textural features 


1. Introduction 

Coronal holes play an important role in geomagnetic storm activity (Tsurutan al., 2 ) and are 

the dominant contributors to space weather disturbances at times of quiet solar activity. The term 
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coronal hole is commonly associated with regions of one dominant magnetic polarity with rapidly 
expanding open magnetic field lines along which solar wind particles escape into interplanetary 
space (Gosling and Pizzo, 1999; Cranmer, 2009). Coronal holes are the source regions of the high 
speed solar wind streams (HSSs). In combination with the Sun’s rotation, they shape the solar 
wind distribution in interplanetary space. In Extreme Ultraviolet (EUV) and X-ray images of the 
Sun, coronal holes are visible as dark areas in the solar corona due to their lower temperature and 
electron density compared to the ambient coronal plasma ( , ). 

The detection of coronal holes purely from their low intensity in solar EUV images is a chal¬ 
lenging task since filament channels also appear as dark coronal features. Filament channels are 
usually interpreted in terms of the weakly twisted flux rope model, having a magnetic field which 
is dominated by the axial component. Dense prominence material is located in the dip of the heli¬ 
cal windings leading to the elongated dark structures observed at a similar dark intensity level as 
coronal holes ( , ; ta, 2011). 

In the past, coronal holes have mostly been identified and tracked by experienced observers. 
There have been recent attempts to automate the process for the identification and detection of 
coronal holes (Barra et al., 2007; Kirk et al., 2009; Krista and Gallagher, 2009; Rotter et al., 2012). 
Those detection techniques are mainly based on differences in intensity compared to the ambient 
corona. 

The extraction of coronal holes is of interest for various applications in astrophysics, e.g. the 
forecast of HSSs. The areas of coronal holes extracted at the solar meridian reveal a correlation 
with the solar wind speed measured ~4 days later in-situ at 1 AU ( , ; 

et al., 2011; Rotter et al., 2012; Rotter et al., 2015). This relation enables us to forecast HSSs based 
on coronal hole observations. Thus, it is important to improve the coronal hole detection method, 
in order to exclude filaments that may erroneously be identified as a coronal hole region in EUV 
images. 

In this work, Solar Dynamic Observatory (SDO) Atmospheric Imaging Assembly (AIA) 19.3 nm 
images ( , ) and SDO/ Helioseismic and Magnetic Imager (HMI) magnetograms 

( al., 2 ) are used in order to analyse the properties of coronal holes and filament chan¬ 

nels during the period 2011 January 1 till 2013 December 31. We first obtain low intensity regions 
from SDO/AIA 19.3 nm images by means of two coronal hole detection techniques: one using 
intensity thresholding ( , ), and a second, called SPoCA ( , ), 

based on fuzzy clustering of intensity values. Comparing these detections with Ha filtergrams from 
Kanzelhohe Observatory (Austria) in which filaments can be clearly identified ( , ), 

we label these regions as either ‘filament’ or ‘coronal hole’. The detected regions were mapped onto 
the corresponding HMI line-of-sight magnetograms. This allows us to calculate a set of attributes 
that takes into account EUV intensity and binary shape information as well as the magnetic field 
structure for each of the regions. 

Our contribution in this paper is two-fold. First, we devise decisive coronal feature attributes. We 
focus on the magnetic flux imbalance, first and second order image statistics and shape measures 
since: (i) coronal holes are expected to be characteristic regions of one dominant magnetic polarity 
whereas filament channels are regions of closed magnetic field lines along a magnetic field reversal 
line; (ii) the physical properties of filament channels are partly reflected in their geometric charac¬ 
teristics. Filament channels are, in contrast to coronal holes, in most cases observable as elongated 
structures that extend in a particular direction. 
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Fig. 1 . Analysis of the detected low intensity regions, (a) AIA 19.3 nm image from 30-Jul-2012 
showing a filament channel; (b) Detected filament channel from EUV image based on intensity 
thresholding; (c) Illustration of the spatial relationship of pixel values contained in the filament 
channel. 


Second, we insert the computed attributes into various supervised classification schemes in order 
to design the most suitable decision rule for a differentiation in real-time. We tested and compared 
through cross-validation four commonly used classifiers, namely Support Vector Machine (SVM), 
Linear SVM, Decision Tree, and Random Forest. To evaluate the performance of different classifiers 
we made use of the Hanssen & Kuipers discriminant, also known as true skill statistics (TSS) 
(Hanssen and Kuipers, 1965; Bloomfield et al., 2012), as a quantitative measure. This paper is 
structured in the following way: In Section 2 the preparation of the labeled data sets is explained. A 
description of the proposed attributes is given in Section 3 and the applied supervised classification 
schemes are presented in Section 4. The results are given in Section 5 and their implications are 
discussed in Section 6. 

2. Data preparation 

Two different coronal hole segmentation techniques (intensity-based thresholding ( , 

; , ), SPoCA ( , )) were applied in order to prepare two data 

sets of low intensity regions during the period 2011 January 1 till 2013 December 31. Based on 
these data we created training sets to design the most suitable decision rule for a differentiation of 
coronal hole and filament channel regions. 

2.1. Coronal hole feature extraction 

2.1.1. Intensity-based thresholding 

An intensity-based thresholding technique was used to detect low intensity regions in SDO/AIA 
19.3 nm images. Based on the intensity distribution in the full-disk EUV images we apply a thresh¬ 
old value of 0.38 x (median intensity value) for areas within ±60° longitude and latitude. Due to the 
optically thin emission from neighboring coronal structures, we used an additional multiplication 
factor of 1.6 for areas outside the ±60° longitude and latitude window. With this refinement, polar 
coronal holes are also well detected although they have a smaller contrast. From this we created 
binary maps, which were post-processed by using morphological methods (erosion, dilation). The 
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effects of these operators on binary maps are to erode (or shrink) the boundaries of regions, remov¬ 
ing all small anomalies (erosion) and to gradually expand the boundaries of regions to fill small 
holes in-between (dilation). A 16 x 16 pixel kernel was used to remove small artefacts and fill small 
gaps in the extracted regions. 

2.1.2. Fuzzy clustering 

The SPoCA-suite ( , ) is a set of segmentation procedures that allows decom¬ 

position of an EUV image into regions of similar intensity, typically active regions, coronal holes, 
and quiet sun. In the framework of the Feature Finding Team project (FFT; ( )), 

the SDO Event Detection System (SDO EDS; ( )) runs the SPoCA-suite to ex¬ 

tract coronal hole information from AIA images in the 19.3 nm passband, and upload the entries 
every four hours to the ffeliophysics Events Knowledgebase ( , ). These entries 

are searchable through the graphical interface iSolSearch, the ontology software package of IDL 
Solarsoft, and the Jffelioviewer visualization tool ( , ). 

The SPoCA extraction method for coronal holes relies on the Fuzzy C-Means algorithm 
(FCM) ( :, ). It clusters the pixel’s intensity values into four classes through the min¬ 

imization of the intra-cluster variance. Such minimization leads to an iterative algorithm, where 
at each iteration every pixel is assigned a membership value between 0 and 1 to each one of the 
class, then the class centers are computed using these memberships. These steps are repeated un¬ 
til convergence of the class centers. The segmented map is obtained by attributing a pixel to the 
class for which it has the maximum final membership value. The coronal hole map corresponds to 
the class whose center has the lowest intensity value. Various preprocessing steps are performed 
before applying FCM algorithm: Images are calibrated using the IDL Solarsoft aia_prep routine, 
and intensities are normalized by their median values. We correct for the limb brightening effect 
by fitting a smooth transition function and inverting it. Finally, a square-root transform is applied 
on the images. Indeed, ( ) showed that for Poisson distributed data a square-root 

transform induces exact asymptotic normality and stabilizes the variance. This is especially useful 
for the extraction of low-intensity regions which are affected by Poisson noise. 

Once the segmented maps are obtained, some post-processing is needed to extract individual 
low intensity regions. Elements with a radius smaller than 6 arcsec are removed and neighbouring 
connected components within a 64 arcsec distance are aggregated. Coronal holes candidates having 
an area smaller than 3000 arcsec 2 are discarded. A study of low intensity feature detected by SPoCA 
during the month of January 2011 revealed that coronal hole candidates detected for more than three 
consecutive days exhibit the expected magnetic properties characteristic of unipolar regions, which 
was not the case for shorter lived regions. Setting up a threshold on lifetime is thus a valuable tool 
to eliminate spurious dark regions, ffence only low-intensity regions which live longer than three 
days are included in the final coronal hole maps used here and in the SDO EDS pipeline (we refer 
to ( ) for more details). 

2.2. Labeled datasets 

The extracted regions were manually labeled as either ‘filament’ or ‘coronal hole’ regions by visu¬ 
ally inspecting Fla filtergrams from Kanzelhohe Observatory (Austria). We selected only regions 
greater than 1900 arcsec 2 for which Kanzelhohe Fla data were available and that were located close 
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to the central meridian (+30° in longitude/latitude). If at the position of an extracted low intensity 
region a filament was clearly observable in the corresponding Ha filtergram (to take into account 
the evolution of filament material we checked a time span of three days), the structure was labeled 
as ‘filament’. Otherwise, the structure was labeled as ‘coronal hole’. This procedure was performed 
for both segmentation methods once per day, resulting in two training sets: one containing 349 
coronal holes and 61 filaments for the intensity-based thresholding method and the second with 
252 coronal holes and 46 filaments for the SpoCA algorithm. In addition we mapped the extracted 
regions onto HMI line-of-sight magnetograms to quantify the magnetic flux characteristics in the 
extracted regions. Note that both segmentation methods have different capabilities for extracting 
low-intensity regions, which explains the dissimilarity in the number of CH candidates detected. In 
future works, we plan on comparing more exhaustively these two methods. 


3. Proposed attributes 

On the prepared training datasets we investigate shape measures, magnetic flux properties and first 
and second order image statistics for the calculation of characteristic attributes. Shape measures 
were used to characterize the shape of the detected low intensity regions. The calculated shape 
measures also allow us to reduce the object pixel configuration contained in a binary map of a 
structure to a single scalar number. With the help of first order and second order image statistics 
(Haralick et al., 1973; Weyn et al., 2000; Ahammer et al., 2008) we focus on the overall pixel value 
distribution and the spatial relation between pixel values contained in AIA EUV images and HMI 
line-of-sight magnetograms. Second order statistics, originally applied on grey scale images, was 
adapted to open value ranges (i.e. including negative data values) in order to calculate textural fea¬ 
tures for intensity and magnetic field configurations. This enables us to characterize the intrinsic 
texture information contained in the extracted structures in AIA 19.3 nm and line-of-sight magne¬ 
tograms via the computation of a set of textural features. A description of the proposed methods is 
outlined in the following subsections. 

3.1. Shape Measures 

While many shape descriptors exist, there exists no generally accepted definition in literature. In 
order to investigate irregular shapes of coronal holes and filament channels in detail, two alternative 
shape measurements were developed. 

Symmetry analysis Shapes are intuitively described with the terms symmetric or asymmetric. We 
measure geometrical symmetry properties with the following technique: After the application 
of discrete geometrical transformations like rotation, reflection and a composition of both, the 
relative overlap in percentage is calculated. 

Direction-dependent shape analysis Shapes can also be presented in form of a function, e.g. the 
average number of neighbours in each direction is representative for the relevant shape informa¬ 
tion of an object. For a detailed description of the used shape measures we refer to 

(2014). 
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3.2. Magnetic Flux Imbalance 

We characterize the magnetic flux imbalance in the detected features by calculating the following 
two attributes from the line-of-sight magnetograms: 

Ri = —, ( 1 ) 

zz_ 


i.e. the ratio of the number of positive (n+) and negative (n_) pixel values, and 


R 2 = 2 


1 

2 


0 + + 10 _| 


( 2 ) 


with 0 + being the total positive magnetic flux in the segmented feature, and ®_ being the total 
negative flux. R 2 quantifies the imbalance between the total positive and the total magnetic flux 
within the extracted regions. In the ideal case, this value would be 1 for coronal holes, as regions 
of a unique polarity (® + = 0 or ®_ = 0), and 0 for filament channels (® + = ®_), as magnetically 
bipolar regions orientated along a polarity inversion line. 


3.3. First Order Image Statistics 


Pixel values contained in the extracted regions were rounded to the nearest integer value. The prob¬ 
ability P(i) for the occurrence of pixels with integer value i is defined as 


P(i) = 


n(i ) 
N ’ 


(3) 


where n(i) is defined as the number of pixels with the actual pixel value i and N is the total number 
of image pixels. The following standard first order attributes were used: 


Mean: /a = ^ i P(i), 

i 

(4) 

Variance: cr 2 = ^ (z - /j) 2 P(i), 

(5) 


i 


Standard Deviation (Contrast): Ci = Vm 2 , 

(6) 

Energy: E l = ^P(i) 2 , 

i 

(7) 

Entropy: S { ^ -P(i) log(P(z)). 

(8) 


Variance cr 2 and contrast Ci are associated with the width of the pixel value distribution. High 
variance or contrast indicate large pixel value differences. Energy E\ is high for unevenly distributed 
pixel values independent of the designated pixel value itself. High entropy S i values reflect the 
degree of information content within the pixel value distribution. Low entropy values correspond 
to coherent pixel values with low information content. High values reflect highly disordered image 
values with a high amount of intrinsic information content. 
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3.4. Second Order Image Statistics 

Second order statistics (Haralick et al., 1973; Weyn et al., 2000; Ahammer et al., 2008) are uti¬ 
lized for the calculation of textural features of AIA and HMI images providing information about 
the spatial arrangement of pixel values. A set of matrices n^jii, j) count how many times a given 
combination of pixel values i and j occur in a particular spatial arrangement. The probabilities of 
co-occurrences of pixel values can then be represented with the so called co-occurrence matrix 
C(i, j). The basic idea how spatial information is contained in low intensity regions is illustrated in 
Figure 1. Textural features allow the description of textural information with a set of scalar numbers. 


3.4.1. Co-Occurrence Matrix 


Analogous to the probability P(i) for the occurrence of pixels with value /', we define P^dii, j ) 
indicating the probability for the co-occurrence of pixel value i and pixel value j in a given distance 
d and direction (p. n^ d (i, j ) is the total number of co-occurrences of pixel value i with pixel value j 
within the distance d and direction (p. The matrix P^dii, j ) is defined as 


j) 


n (K d(i, j) 
N 


(9) 


where N is the total number of object pixels. Since there is no consensus about the choice of angles 
(p and distance d, the calculations have been done with a distance of one pixel and 8 directions 
considering only the nearest neighbourhood of each pixel. The used set of angles Q = { <p 0 ,..., (pi) 
is given by 


(pi — l • 45 , l — 0,... ,7 , (pi G D. 


( 10 ) 


This set of angles provides 8 matrices P^d(i,j), one matrix for each of the directions. The co¬ 
occurrence matrix C(i, j ) is finally defined as 

C(J, j ) = avg (P (l) Ji, j)). (11) 

The entry (z, j) in the co-occurrence matrix C(z, j) therefore indicates the probability for the co¬ 
occurrence of pixel value i with pixel value j within the nearest neighbourhood. In Figure 2 the 
corresponding co-occurrence matrix for a filament channel, shown in Figure 1, is visualized. For 
the implementation it is necessary to differentiate between the actual pixel value and the position in 
the co-occurrence matrix. 

Important properties of the co-occurrence matrix C(z, j) are represented in Figure 2: (i) C(z, j) is 
always a symmetric matrix since the number of occurrences for pixel value i together with pixel 
value j will be always the same as vice versa, (ii) C(z, j) provides a probability density function 
(PDF), to find probabilities of joint relationships between a pair of pixel values i and j. Thus, the 
sum over all entries must be equal to 1. 


3.4.2. Textural Features 

A set of textural features (Appendix A) can be calculated from the co-occurrence matrix C(z, j). For 
illustrative propose, we discuss two textural features (energy H { and contrast H 2 ) in detail: 

w i = ZZ C(i ’ i)2 ’ (12) 

i j 
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Fig. 2. Representation of the co-occurrence matrix C(i, j ) calculated from the filament channel shown in 
Figure 1. The colorbar of the co-occurrence matrix entries C{i,j) represent the probability for neighbouring 
occurrences of pixel value i with pixel value j. 

= C(i,j). (13) 

i j 

Energy H\ is a measure of homogeneity of the pixel pair distribution. It is high for a small number 
of different pixel pair combinations and low for a large number of different combinations. H\ is at 
its maximum if the image contains only a single pixel value, independent of the pixel value itself. 
Thus, the energy is mainly a measure of the sharpness of the peak in the co-occurrence matrix. It 
is high for a sharp peak, indicating only a few pixel value combinations, and low for a broad peak, 
indicating a variety of pixel value combinations. 

Due to the quadratic influence of the differences in the contrast H 2 , the value is high for large 
differences between pixel value i and pixel value j. Hence, mainly the position of the peak in the 
co-occurrence matrix is important (see Figure 2). The position indicates small or large differences 
of the pixel values. The contrast feature is a measure of the amount of local variations present in the 
image. It is high for large differences of pixel pair values and low for small differences. All textural 
features that were used in this analysis are listed in Appendix A. 

4. Supervised classification problem 

The computed attributes were used as input for data mining in order to design the most suitable 
decision rule for a separation between coronal holes and filaments. Supervised classification is 
commonly applied when a large set of attributes is available and the interpretation of the obtained 
information is complicated. In supervised classification problems, an object is observed and we 
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aim to classify it into one of two classes (0 or 1, here ‘filament’ or ‘coronal hole’). To make this 
decision we have access to measurements of various properties of the object, called features. Given 
these feature vectors x 6 R d , we aim to find a decision (or classification) rule, that is, a mapping 
c : R d —> {0,1}, where c(x) indicates the decision when feature vector x is observed. Out of the 
numerous collection of possible decision rules, we aim to identify those performing best. Four 
widely used algorithms with proven theoretical properties (Section 4.1) are evaluated via a common 
protocol (Section 4.2). Section 5 describes the performance measures used to compare them. 

The set of attributes described in Section 3.1-3.4 together with the labelling of the maps as ‘fila¬ 
ment’ or ‘coronal hole’ provide us with a labeled dataset of feature values. We consider two situa¬ 
tions. First, the whole set of attributes, computed with both SDO/AIA and SDO/HMI from SPoCA, 
is used as feature values with the aim to reach the best classification. Second, we create new clas¬ 
sifiers trained on a sub-set of attributes using only AIA information on SPoCA maps. The aim is 
twofold: (i) obtain a set of rules that are easy to implement into the SDO EDS pipeline, (ii) measure 
the corresponding improvement in performance for the SPoCA algorithm. 

4.1. Supervised classification algorithms 

We tested four classifiers using the scikit-learn library ( , ): 

Linear Support Vector Machine The key idea of Linear Support Vector Machine (Linear SVM) 
is to find hyperplanes that separate the data as much as possible, that is, with a large margin. 
SVM optimizes a trade-off between maximizing the margin of separability between classes and 
minimizing margin errors. It provides a convex approximation to the combinatorial problem of 
minimizing classification errors. The practical implementation is formulated as the minimization 
of a penalized loss function. We used the LIBLINEAR ( , ) implementation of 

Linear SVM. 

Support Vector Machine A second key idea of SVM as presented by Vapnick in his original for¬ 
mulation ( , ) is to map the feature vectors in a nonlinear way to a high (possibly infi¬ 

nite) dimensional space and then utilize linear classifiers in this new space. This is done through 
the use of a kernel function. In our case, we tried several kernel functions: gaussian, sigmoid, 
polynomial and linear from the LIBSVM library ( , ). 

Decision Tree Decision trees produce a set of if-then-else decision rules that are selected on the ba¬ 
sis of the expected reduction in entropy that would result from sorting a particular attribute. Trees 
are usually grown to their maximum size before a pruning step is applied to reduce over-fitting. 
Various decision tree algorithms have been proposed ( , ; , ). Lor 

this work we used decision trees using the Gini and entropy criteria and various depths and split¬ 
ting rules 1 . 

Random Forest We tested an ensemble classifier method 2 . More precisely, we used Random 
Lorest, where a set of decision trees is created by introducing some randomness in the con¬ 
struction of the decision tree. The prediction of the ensemble is given as the averaged prediction 
of the individual decision tree ( , ). 

1 http://scikit-learn.org/stable/modules/tree.html 

2 http://scikit-learn.org/stable/modules/ensemble.html 
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4.2. Training and evaluation protocol 

The supervised classification in machine learning is always performed in two steps. 1.) During the 
training phase, a model is estimated from the data. 2.) In the validation (or test-) phase the trained 
model is applied on other data and model properties (such as error classification rate) are estimated. 
This is done in practice by splitting the initial data set into a ‘training dataset’ and a ‘validation 
dataset’. 

However, in the present study we want to compare the performance of several classifiers. In this 
case, performing the splitting only once is not enough, as the observed difference between two 
classifiers may depend on the chosen training and test samples. To avoid this, we need to repeat 
the comparison over randomly selected partitions of the data and report the average performance 
(see Chap. 5 in ( ) for a discussion on error estimation in classification 

problem). 

We therefore perform 100 iterations of the following protocol. We do a stratified shuffle-split of 
the original dataset into a 75% development set and 25% evaluation set. This means that we shuffle 
the original dataset, then split it into two 75% / 25% subsets (shuffle-split) where each subset has 
approximately the same class distributions as the full dataset (stratification). The development set 
in each iteration is used to train and evaluate each hyper-parameter combination for each algorithm. 
To choose the best hyper-parameter combination we use stratified 5-fold cross-validation. A A-f'old 
cross-validation means that the development set is further split into k-folds. Each combination of 
(k - 1) folds is used for training and the remaining fold serves for testing, for a total of k train/test 
splits. The use of a stratified 5-fold cross-validation means that each fold has approximately the 
same class distributions as the original dataset. 

Once the optimum combination of hyper-parameters is found, it is used to train a classifier on the 
75% development set, and is evaluated on the 25% evaluation set. This final hold-out evaluation set 
is necessary to accurately estimate real-world performance because the cross-validation can over-fit 
for the particular split. The performance is measured by computing a skill score for each of the 100 
iterations, see Section 5. This allows us to quantify an average performance, but also to evaluate the 
variance in performance results across different runs. The analysis code can be accessed online 3 . 


5. Results 

To evaluate the performance of the four classifiers we computed for each of the 100 runs the ‘confu¬ 
sion matrix’. A confusion matrix, see Table 1, contains the elements TP (true positive, coronal hole 
predicted and observed), FP (false positive, coronal hole predicted and filament channel observed), 
FN (false negative, filament channel predicted and coronal hole observed) and TN (true negative, 
filament channel predicted and filament channel observed). 

There exists a variety of skill scores (Woodc , ) built by combination of the confusion 

matrix elements. One of the most frequently used is the Hanssen & Kuipers discriminant also known 
as true skill statistics (TSS) (Hanssen and Kuipers, 1965; Bloomfield et al., 2012). The TSS is 


3 https://github.com/rubendv/ch_filament_classification 
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Predicted 

Observed: Coronal Hole (CH) 

Filament Channel (FC) 

CH 

True Positive (TP) 

False Positive (FP) 

FC 

False Negative (FN) 

True Negative (TN) 


Table 1. Coronal hole and filament channel classification contingency table (confusion matrix). 


defined as the proportion of correctly predicted coronal holes among all coronal holes (TPR) minus 
the proportion of filaments that were classified as coronal holes among all filaments (FPR): 4 


TSS = TPR - FPR = 


TP 

TP+FN 


FP 

FP+TN' 


(14) 


The TSS is defined in the range [-1,1], A TSS of 0 indicates that the algorithm cannot distinguish 
between coronal holes and filament channels. A perfect classifier would have the value 1 or -1 
(inverse classification), respectively. In ( ) the TSS was proposed to be the 

standard skill score for comparing the performances of flare forecasts with differing flare/no-flare 
sample ratios and is also appropriate for this study. When applying the developed algorithms in the 
form of an automated coronal hole detection tool we are mainly interested in both, the proportion 
of correctly predicted coronal holes (TPR) as well as the proportion of falsely predicted coronal 
holes (FPR). Hence, the combination of TPR and FPR in form of the TSS provides an important 
and intuitive indicator for the performance of the developed framework. Furthermore, the TSS uses 
all elements of the confusion matrix and is the only known skill score which is unbiased by the 
proportion of coronal holes and filament channels in the data sets. 

Figure 3 (resp. 4) shows for the intensity based threshold method (resp. for the SPoCA method) 
the proportion of correctly predicted coronal holes (TPR) (that is, the efficiency or sensitivity) versus 
the proportion of falsely predicted coronal holes (FPR). The performance for each of the 100 runs is 
visually represented in those figures by the points (FPR, TPR). A perfect classification would have 
only points at the position (0,1), a perfect inverse classification points at (1,0) and ‘no performance’ 
is represented by the diagonal line (FPR = TPR). 

Figure 3 shows the performance of the four classifiers using all attributes (from AIA and magne¬ 
tograms) computed with the intensity-based threshold segmentation maps. All classifiers produce 
a high TPR, but only the two SVMs achieve at the same time a consistently low FPR, with Linear 
SVM performing slightly better than SVM. In Figure 4 the results for the SPoCA detection of 
SVM and Linear SVM are presented. Both classifiers achieve high TPR and low FPR. A compari¬ 
son of the performance with and without magnetogram information indicates that the inclusion of 
magnetogram-based attributes in addition to the AIA-based attributes systematically improves both 
TPR and FPR measures across all four classifiers: Figure 4 compares the density of TPR versus 
FPR obtained for both SVM classifiers when adding or not the magnetogram information, Table 2 

4 For the reader who is used to the terms ‘recall’ and ‘precision’ we note that recall is equivalent to the TPR, 
and precision denotes the fraction of correctly predicted coronal holes among all predicted coronal holes. The 
TSS is given by TP/(TP+FN) + TN/(TN+FP) - 1, the sum of recalls for coronal holes and filament channels 
minus a scaling factor of 1. 
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gives the median TSS and the standard deviation for all tested classifiers, and Figure 5 provides a 
graphical comparison of the actual TSS values in the form of box plots. 




True Skill Statistics (TSS) 


Datasets 

SVM 

Linear SVM 

Random Forest 

Decision Tree 

Intensity Threshold (all attributes) 
SPoCA (all attributes) 

SPoCA (only attributes from AIA) 

0.90 ± 0.07 

0.89 ± 0.07 
0.82 ± 0.08 

0.94 ± 0.05 

0.90 ± 0.06 
0.83 ± 0.07 

0.87 ±0.10 

0.77 ±0.11 
0.71 ±0.11 

0.80 ±0.15 
0.68 ±0.18 
0.63 ±0.21 


Table 2. Median TSS and standard deviation for all four classifiers and the three different datasets. 

For both segmentation techniques best results are achieved with linear SVM showing the highest 
median TSS and smallest standard deviation. Similar results can be obtained with SVM, although 
with slightly lower values of the median TSS. In contrast, Decision Tree and Random Forest exhibit 
significantly lower median TSS and higher standard deviation. For the dataset extracted with the 
SPoCA algorithm we also considered the case where attributes are obtained from AIA images 
only. The corresponding performance is, as expected, lower than when magnetogram information 
is included. 

Another comparison of the investigated classifiers is presented by the Receiver Operator 
Characteristic (ROC) curves. ROC curves are commonly used to compare and evaluate the per¬ 
formance of different classifiers. They characterize how the number of correctly classified coronal 
holes (TPR) varies with the number of incorrectly classified filament channels (FPR). This gives 
an illustration of the tradeoff between the completeness of coronal holes and the contamination 
with filaments for each classifier. Hence, the ROC curves convey the information of which classi¬ 
fier should be applied when either a small, focused sample of coronal holes or a larger sample of 
coronal holes is requested. 

In Figure 6, the ROC curves for SVM and linear SVM calculated from the combined results of 
each of the 100 iterations are shown. Again, the FPR is plotted on the x-axis and the TPR on the 
y-axis. The best results are represented by curves close to the upper left comer with a possible large 
area under the curve (AUC). The AUC is usually used in order to roughly summarize the obtained 
ROC curves in a single quantity. The best performances, indicated by a high completeness together 
with a low contamination, are achieved with the Linear SVM followed by the SVM using all com¬ 
puted attributes. The ROC analysis confirms that the exclusion of the magnetogram information 
negatively affects the performances of both classifiers by decreasing the TPR and increasing the 
FPR. 

6. Discussion and Conclusion 

The most popular methods for the extraction of coronal hole regions are based on the intensity 
information in EUV images of the Sun. However, other features, notably filament channels, ap¬ 
pear similarly dark as coronal holes. It is a general problem of intensity-value-based techniques 
to distinguish between those features. We found in our carefully prepared datasets using two dif¬ 
ferent segmentation algorithms (intensity-based thresholding and SPOCA) that about 15% of the 
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FPR FPR 




Fig. 3. Density plot of correctly predicted coronal holes (TPR) versus erroneously predicted coronal 
holes (FPR) obtained from intensity-based coronal hole region segmentation using all attributes. 
The red crosses denote the individual results for the TPR and FPR for each of the 100 iterations, 
according to the training and evaluation protocol described in Section 4.2. (a) Decision Tree; (b) 
Linear SVM; (c) Random Forest; (d) SVM. 


coronal hole candidates are actually filament channels. We therefore developed new algorithms and 
techniques to tackle this shortcoming. 

We investigated a new approach using a combination of segmentation techniques and ma¬ 
chine learning algorithms. The proposed attributes were used as input for supervised classifica¬ 
tion schemes. Four commonly applied classifiers, including SVM, Linear SVM, Decision Tree and 
Random Forest were tested for two different segmentation techniques and evaluated based on true 
skill statistics. The results reveal that all classifiers provide good results in general, with the linear 
SVM and SVM classifier providing the best performances (TSS ~ 0.90). We also found that the two 
different image segmentation techniques provide very similar results indicating the broad applica¬ 
bility of the presented method. In addition to EUV images we used the magnetic field information 
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Fig. 4. Density plot of correctly predicted coronal holes (TPR) versus erroneously predicted coronal 
holes (FPR) obtained with SPoCA algorithm. The red crosses denote the individual results for the 
TPR and FPR for each of the 100 iterations, according to the training and evaluation protocol 
described in Section 4.2. (a) SVM (all attributes); (b) Linear SVM (all attributes); (c) SVM (AIA 
attributes only); (d) Linear SVM (AIA attributes only). 


from HMI line-of-sight magnetograms. Adding the information and extracted attributes from the 
magnetic field data improves the performance across all four classifiers for the SPoCA detection 
(c.f. Table 2). 

Besides standard statistical attributes we investigated the benefits from textural features, origi¬ 
nally introduced by ( ), to analyze the texture information of coronal holes and 

filament channels in high spatial resolution AIA 19.3 nm and line-of-sight magnetograms. For this 
purpose we adapted the original method introduced by Haralick to open value ranges in order to 
calculate a set of textural features for intensity and magnetic field configurations. This enabled us 
to describe the intrinsic spatial arrangement contained within coronal hole and filament channel 
regions in AIA 19.3 nm and line-of-sight magnetograms with a set of textural features. Textural 
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Fig. 5. Comparison between the computed true skill statistics (TSS) for different classifiers using 
different segmentation techniques and sets of attributes. On each box, the central mark is the median 
(50th percentile), the edges of the box are the 25th and 75th percentiles, the whiskers show the 
±2.1cr range covering 99.3% of the data and outliers are plotted individually as red crosses, (a) 
Intensity Threshold (all attributes) (b) SPoCA (all attributes); (c) SPoCA (only attributes from AIA). 
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FPR 


Linear SVM (SPoCA, all attributes), AUC = 0.986 

- Linear SVM (SPoCA, attributes from AIA), AUC = 0.970 

- SVM (SPoCA, all attributes), AUC = 0.980 

- SVM (SPoCA, attributes from AIA), AUC = 0.971 

Linear SVM (Int. Threshold, all attributes), AUC = 0.963 
- SVM (Int. Threshold, all attributes), AUC = 0.960 


Fig. 6. ROC curves for the SVM and Linear SVM classifier together with the calculated area under 
the curve (AUC). 


features provide information about the structural arrangement of pixel values in contrast to first 
order statistics, which focus on the overall pixel value distribution contained in an image. In this 
respect we note that for the detailed analysis of textural features one has to take into account the 
influence of noise. Digital images are discrete representations of objects or scenes, unavoidably, 
they have limited spatial resolution and they are contaminated with noise ( , ). 

Hence, the image quality is a limiting factor for using such ki nd of methods and low signal to noise 
ratios may hinder the application of textural information. 

The algorithms are not costly with respect to computing-time. The segmentation algorithms as 
well as the calculation of attributes and classifiers are calculated within a few minutes. This is 
favorable for real-time application. We therefore aim to implement the described algorithms in 
currently available forecasting tools for solar wind high speed streams 5 . We also plan on improving 
the SPoCA-CH module within the SDO EDS pipeline using Linear SVM classification on AIA 
attributes. 

In summary, this study demonstrates how a machine learning approach may help improve upon 
an unsupervised feature extraction method 6 . The findings suggest that the proposed attributes and 
classifiers may actually be applicable for a wide range of imaging data. We stress that the inclusion 
of magnetic field data improved the performance of the coronal hole and filament discrimination in 
EUV images. All tested classifiers provide good results for the TSS and we therefore expect that the 
developed detection tool has the potential to reduce coronal hole classification error rate for both 
segmentation algorithms. 
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Appendix A: Textural Features 

A.l. Notation 

We suggest to use a set of textural features which can be calculated from the co-occurrence matrix. 
The following equations define these features. 

C(i, j) (i, j)-th entry in the normalized co-occurrence matrix. 
p x (i) i-th entry is obtained by summing the rows of C(i, j). 

p y (j ) j-th entry is obtained by summing the columns of C(i, j). 

p x -y(k) k-th entry of p x - y (k) corresponding to the sum over all entries of C(i, j ) with 
absolute pixel value difference of i and j equal to k. 
p x+y {k ) k-th entry of p x+y (k) corresponding to the sum over all entries of C(i, j ) where 
the addition of i and j is equal to k. 
p x ,p y means of p x and p y . 

<t x , (T y standard deviations of p x and p y . 

N g Number of distinct pixel values. 

'Lh'Lj Convention indicating • 

HX, HY Entropies of p x and p y 


Px® = 2 J C ®j )’ . PyU) = ^C{i,j), 

j i 

(A.l) 

p x+y (k ) = ^ J) > i + j = k, k = 2,3,..., 2 N g , 

i j 

(A. 2) 

p x -y(k) = ^ X C<yi ’ 7 ), \i - j\ = k ’ k = 0,1,..., Ng - 1, 

*' j 

(A.3) 

HX = - ^ p x (i) log (p x (i)) , HY = - py(i) log (p y (i)) 

i i 

(A. 4) 

HXYl =~YjYj C(U j) lo § (Px®PyU)), 
i j 

(A.5) 

HXY2 = - J] J] Px® PyU) log (Px® PyU)), 

(A. 6) 


i j 


A. 2. Textural Features 

Based on this notation, the following textural features can be calculated: 

Energy: 

Hi (A.7) 

i j 
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Contrast: 


i j 

(A.8) 

Correlation: 


H O' _ Px )0 “ Py) 

i j tT x f - r y 

(A.9) 

Sum of Squares - Variance: 


H 4 = fifCd, j), 

i j 

(A. 10) 

Homogeneity: 


tt _ V V 

5 ZjZj 1 + (i-j) 2 ’ 

1 J J 

(A. 11) 

Sum Average: 


2N g 


Px+y(l), 

;=2 

(A. 12) 

Sum Variance: 


2N g 


= d - H s ) 2 p x+y (i), 

i= 2 

(A. 13) 

Sum Entropy: 


2N g 


Hz = - ^ p x+ y(i ) log (p x+y (/)), 
i=2 

(A. 14) 

Entropy: 


^9 = -XZ C(/ ^' )1 ° g(C(/ ^' )) ’ 

I j 

(A. 15) 

Difference Variance: 


//to = var (/t A _ y ) , 

(A. 16) 

Difference Entropy: 


N s ~ 1 


#11 = - J] /fv-v(') log(p A -v(/)), 

j =0 

(A. 17) 

Information Measures of Correlation (I): 

H 9 - HXY1 
max(HX,HY) 

(A. 18) 

Information Measures of Correlation (II): 

H n = V1 - exp(-2 (HXY2 - H 9 ) ). 

(A. 19) 
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